home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / ChiSquare.c < prev    next >
Text File  |  1995-05-27  |  1KB  |  58 lines
  1. /*
  2. ChiSquare.c
  3. Copyright 1990-1992 © Denis G. Pelli
  4. PChiSquare() calculates the significance level at which a fit may be rejected.
  5.  
  6. Also see: Binomial.c, Exponential.c, Normal.c, Uniform.c
  7.  
  8. HISTORY:
  9. 1/7/84    dgp    wrote it in FORTRAN
  10. 8/9/89    dgp    translated it to C
  11. 4/5/90    dgp    tested it against tables in Abramowitz in Stegun. Works fine.
  12. 4/8/90    dgp    renamed file ChiSquare.c from PChiSquare.c
  13. 7/30/91    dgp    now use NAN defined in VideoToolbox.h
  14. 1/12/92    dgp    Renamed NormalPDF() to NormalPdf().
  15. */
  16. #include "VideoToolbox.h"
  17.  
  18. double PChiSquare (double chiSquare,int n)
  19. /*
  20. Returns one minus the cumulative probability of the
  21. Chi-Square distribution for value chisq with n degrees of freedom.
  22. The formula is from Abramowitz and Stegun, Eqs. 26.4.4 and 26.4.5, pg 941.
  23. The formula is exact.
  24. Range: chisq>=0.0, n>=0,
  25. Returns zero if n is zero.
  26.  
  27. M. Abramowitz and I.A. Stegun (1964) Handbook of mathematical functions. Dover.
  28. */
  29. {
  30.     double x,a,P;
  31.     int i,ii;
  32.     
  33.     if(chiSquare<0.0 || n<0) return NAN;    /* return NAN to signal error */
  34.     x = sqrt(chiSquare);
  35.     if(n%2 != 0) {
  36.         /* n is odd */
  37.         P = 2.0*Normal(-x);
  38.         a = 2.0*NormalPdf(x)*x;
  39.         ii = (n-1)/2;
  40.         for(i=1;i<=ii;i++){
  41.             P += a;
  42.             a *= chiSquare/(2*i+1);
  43.         }
  44.     }
  45.     else {
  46.         /* n is even */
  47.         P = 0.0;
  48.         a = sqrt(2.0*PI)*NormalPdf(x);
  49.         ii = n/2;
  50.         for( i=1;i<=ii;i++){
  51.             P += a;
  52.             a *= chiSquare/(2*i);
  53.         }
  54.     }
  55.     return P;
  56. }
  57.  
  58.